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This paper presents the overall mathematical model and results from pseudo linear 
recursive estimators of attitude and rate for a spinning spacecraft. The measurements 
considered are vector measurements obtained by sun-sensors, fixed head star trackers, 
horizon sensors, and three axis magnetometers. Two filters are proposed for estimating 
the attitude as well as the angular rate vector. One filter, called the q-Filter , yields the 
attitude estimate as a quaternion estimate, and the other filter, called the D-Filter , 
yields the estimated direction cosine matrix. Because the spacecraft is gyro-less, Euler’s 
equation of angular motion of rigid bodies is used to enable the estimation of the 
angular velocity. A simpler Markov model is suggested as a replacement for Euler’s 
equation in the case where the vector measurements are obtained at high rates relative 
to the spacecraft angular rate. The performance of the two filters is examined using 
simulated data. 


I. Introduction 

The advantages of inducing a constant spin rate on a spacecraft are well known. A variety of science 
missions have used this technique as a relatively low cost method for conducting science. Starting in the later 
part of the 1970s, NASA pushed toward building spacecraft using 3-axis control as opposed to the single-axis 
control mentioned above. Considerable effort was expended toward sensor and control system development 
as well as the development of ground systems to independently process the data. As a result, spinning 
spacecraft development and their resulting ground system development stagnated. In the 1990s, shrinking 
budgets made spinning spacecraft a more attractive option for science though most NASA missions 
continued to be three-axis stabilized. The attitude requirements for spinning spacecraft have become more 
stringent and the ground systems must be enhanced in order to provide the necessary attitude estimation 
accuracy. 

The current ground attitude determination system at NASA Goddard Space Flight Center consists of a 
least squares estimator that assumes constant dynamics. Since spinning spacecraft (SC) have no gyroscopes 
for measuring attitude rate, any new estimator would need to rely on the spacecraft dynamics equations. One 
estimation technique that utilized the SC dynamics and has been used successfully on 3-axis spacecraft 
ground systems is the pseudo-linear algorithm. Therefore in this paper, too, SC dynamics as well as spinning 
spacecraft sensor models are developed to work with the pseudo-linear attitude and rate estimation algorithm. 
In the rest of this section we present the definition of the coordinate systems used in the filter developments 
and explain how the important angles are extracted from the estimated variables. In the following section we 
develop the sensor measurement mathematical models, then, in Section III, we present two filters for attitude 
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determination. The filter algorithms are summarized in Section IV, the test that was done is presented in 
Section V and the conclusions are summarized in the last section. 

1.1 Definition of Coordinate Systems 

The transformation from the inertial coordinate system to the spacecraft (SC) body system is expressed by 
a sequence of three Euler angles a , (3 , and y in the order 3-2-3. This sequence is illustrated in Fig. 1 



Fig. 1: 3-2 -3 transformation by a, (3, y from inertial to body coordinates. 


where the superscript / denotes the inertial coordinates. The first rotation by a about the Zj axis yields a 
coordinate system denoted by 1. The next rotation about Yj by the angle [3 results in the intermediate 
system denoted by 2. Note that Z 2 is identical to the body axis, Z b , which is the spin-axis. One more 
rotation by y about this spin axis yields the body coordinate system denoted by b. Actually y is the spin- 

angle, which is also known as phase-angle. For this sequence of Euler angles the direction cosine matrix 
(DCM) which transforms vectors from inertial to body coordinates is: 
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SC Z b - axis. From the geometry it is clear that the three components of the unit vector to the sun expressed 
in body coordinates are as follows* 


z b 



Fig. 2 : Sun vector geometry in body coordinates 


s bx = scp-cS 

s by =scp-sS (3. a) 

s b2 = cep 

where S bx , S by and s bz are the three components of the unit vector to the sun, l sb expressed in body 

coordinates. Note that $ is a fixed angle and only (p is changing. 

If the sun sensor coordinates are not identical to the body coordinates then the measured vector is given in 
the sun sensor coordinates and it has to be transformed to the SC body coordinates as follows. 


Although the sun sensor may measure functions of the angles and not the angles themselves, for simplicity we assume that the angles 
are the measured quantities. It is indeed easy to express the sine and cosine functions of the respective angles in term of the measured 
quantities. 
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where I is the SC inertia tensor, CO is the angular velocity vector, h is the angular momentum of the 
momentum wheels, and T is the external torque acting on the SC. The symbol [a x ] denotes the cross 
product matrix of the general vector a . Attitude is represented by the attitude quaternion whose dynamics 
equation is [Ref. 1, p. 512] 


where 



( 8 ) 


( 9 ) 


In the q-filter we augment Eqs. (7) and (8) to form the following dynamics equation, which includes the 
noise terms and W q . 
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The unbiased white-noise vector accounts for the inaccuracies in the modeling of the SC angular 
dynamics, and W q is an unbiased white-noise vector that accounts for modeling errors in the quaternion 
dynamics. 

When the measurements come in at a relatively high frequency we may be able to replace the SC angular 
dynamics model in Eq. (10. a) with a simpler Markov model 7 . Consequently Eq. (lO.a) is replaced by the 
model 
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where [x] is a diagonal matrix whose elements are the inverse of suitable time constants. 

III.l .2 The Measurement Model 

Let I*j denote the column vector measured by the sensor and expressed in the reference inertial 
coordinate system (this vector is taken from the almanac). The relationship between this column vector and 
b j , which is that vector resolved in body coordinates, is expressed by 

b-Dlr (11) 


The measured bj vector, which we denote by b jm , contains an additive zero-mean white-noise vector, 
v jib , thus 
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jm 
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Using the expression for given in Eq. (2.b), and after some algebraic manipulation, the last equation can 
be written as 


where 


h jM>- 
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(13 a) 


(13. b) 


It should be noted that H is not unique; that is, D^(q)r. can be written in the shape H.(r.,q)q where 
H.(r.,q) can take several forms. However, the form of H.(r.,q) shown in Eq. (13.b) is the most symmetrical 

one and it is exactly half the sensitivity matrix, §[D^ (q)r ] 6q , obtained when developing the corresponding 
Extended Kalman Filter measurement equation. Finally, Eq. (13. a) can be written also as 


b ,-[° H ,M>] 



(13.C) 


The last equation is the measurement equation associated with any vector measurement. 

III. 1.3 Vector Propagation 

As mentioned before, it was found in the past 6 that in non-spinning SC a similar filter performed well 
only when the different vector measurements were taken at short time intervals from one another. Therefore 
we propose to propagate forward each measurement after it has been used, and use the propagated vector 
together with the newly acquired vector measurement. The projection is done as follows. Suppose that at the 

time when b 1 is measured, the quaternion of the transformation from inertial to body coordinates is q^, 
and at the time when b 2 is measured, that quaternion is q Ib2 , then q Ib2 satisfies the following relation 


q, b2 =q*,®5q 


(14. a) 


where 5q is the quaternion that transforms from the body coordinates at the time b : is measured to the 
body coordinates at the time b 2 is measured. From Eq. (14. a) we obtain 


Sq = q, b r‘®q Ib2 


The corresponding DCM, which we denote by 5D' 2 , is computed using 
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hence the projected b, , denoted by b!’ is computed as 
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(14. d) 


b>5D;b, 


Both bj 5 and b 7 are used in the filter update at the time b 2 is measured. If bj is a sun sensor or a star 
tracker measurement, then rf , the vector corresponding to bj 3 , does not change between the two 
measurement time points; that is, r/ 5 = r, . However, when the measurement is of the magnetic field or of 

the nadir, does change, but the new vector is easily computable independently of the attitude. 

When this approach is used, situations may arise where the attitude is unobservable even though attitude 
is completely defined from two pairs of vectors [see Ref 6 for such a case]. To avoid such cases it is 

proposed to use the technique that was used successfully in 6 . Generate a pseudo-measurement pair, r 3 and 
b 3 where 


and 


= r , P X r . 

(15.a) 
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Using Eq. (12) the noise vector, v 3 , associated with the pseudo-measurement b 3 is given by 


v 3 =s (SD 2 r 1 )xv 2 +v l x(5D‘r 2 ) + v l xv 2 (16) 

Note that v 3 is correlated with both v ] and v 7 . Since as well as V 7 consist of components that are 
zero mean sequences which are uncorrelated with themselves at different time points and with each other at 
all time points, it is not too difficult to compute the covariance matrix of R 3 ; however, since the suitable 
measurement covariance matrix is normally determined experimentally by tuning 3 , the exact computation 
of R 3 is, in practice, unnecessary. The corresponding measurement model when using vector projection is 
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where the subscript m denotes measured quantity. The covariance matrix that corresponds to SD^Vj is 
SD.R^D^) 1 where R, is the covariance matrix corresponding to V ] . Because of the relationship, 

expressed in Eq. (16. a), between V 3 and Vj , and between V 3 and v 2 , it is obvious that the covariance 

matrix of the noise vector in Eq. (17) has off-diagonal elements. The simplest way to avoid the computation 
of these elements is to disregard these elements and compensate for their elimination by proper tuning. 

III.2 The D-Filter 
III.2.1 The Dynamics Model 

We start with the well-known first order DCM differential equation 9 

D = -[cdx]D (18) 

where [co x] is the cross product matrix of CO defined as follows. 
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Define the columns of D as 3 vectors; that is, 


D = [d, |d 2 |d 3 ] 


then Eq. ( 1 8) can be written as 


It is easy to see that 
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Define the matrix ® as follows 
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where, obviously, 0 € 5f 9x3 . Also define the vector <£ as follows 



then Eqs. (22) can be written in one vector differential equation as 
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In order to use the last equation as a dynamics equation in a KF, we need to add to it zero-mean white-noise 
to account for modeling inaccuracies. This results in the following dynamics equation. 
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where, as before, W ffl accounts for the inaccuracies in the modeling of the SC angular dynamics and 

accounts for modeling errors of the DCM dynamics. As with the q-filter, here too, when the measurements 
are obtained at a relatively high frequency, we can use the simpler dynamics equation. In such case the 
dynamics equation reduces to 
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III.2.2 The Measurement Model 

Equation (12) can be written as 


b ;. m = [ d 1 r i l +d 2 r j 2 +d 3 r u] + v i.b (29) 

which can be written also as 

b = R./ + v (30) 

j.m j j,b v } 

where 

R j = [ r i. r 3 fj, 1 , (31) 

and (£ is as defined in Eq. (25). Note that usually V- is taken from an almanac and is not measured; 
therefore, it is very accurate. The measurement model is 
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CO 
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III.2.3 Vector Projection 

As with the q-fllter, here too we propose vector projection. Here, however, SD^ is obtained as follows 


D 1 =5D 1 D I 

b2 2 b! 


(23. a) 


where D^, is the transformation from inertial to body when the measurement is of b, and is the 
transformation from inertial to body when the measurement is of b 7 . From Eq. (33. a) we obtain 


5d;=(d;,) t d; 2 (33.b) 

The handling of the vector measurements in this filter is identical to that of the q-filter. The measurement 
model is different though. In the D-filter it is as follows. 
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IV. SUMMARY 


IV.l Using the q-Filter 

For the dynamics equation use either 


CO 


T'[(Ie>+h)x] o’ 

CO 

+ 

’r'(T-h)' 

+ 

W 

w 

A. 


iQ o_ 

_q. 


0 


W 

L q J 


CO 



'-M 

o' 

co~ 


'o' 

+ 

W 

w 

_q. 


.tQ 

0 _ 

.qj 


0 


w 

L q J 


11 

American Institute of Aeronautics and Astronautics 


(lO.a) 


(lO.b) 



and for the measurement equation use 
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V. TEST RESULTS 

V.l The q-Filter Results 

In this simulation we considered a SC spinning at a rate of 120 deg/sec. The initial Euler angles 
describing the attitude were <2 =- 86 . 4 , (i ~ 33.3 ,y = 8.2 degrees. The simulated measurements 

included Sun Sensor (SS) measurements, Three Axis Magnetometer (TAM) measurements and Sun Pulse 
time. The latter was used to obtain a rough initial spin rate estimate. For better convergence, the initial SS and 
TAM measurement were used in the TRIAD algorithm 10 to obtain a rough initial attitude estimate. The 
attitude and rate errors mean and standard deviation are presented, respectively in Tables I and II. Plots of the 
attitude and rate errors are shown in Figs. 3 and 4. 


Table I: Attitude error mean and standard deviation using the q-filter (deg). 



Component 

Mean error 

Standard Deviation 

ATTITUDE 

ERRORS 

X 

-0.068204 

0.036345 

y 

-0.003806 

0.064319 

z 

0.076284 

0.002280 


Table II: Rate error mean and standard deviation using the q-filter (deg/sec). 



Component 

Mean error 

Standard Deviation 

RATE 

ERRORS 

X 

-0.244370 

0.076580 

y 

-0.007258 

0.011394 

z 

0.000006 

0.000113 


Spinner Attitude Errors (deg) 



Fig. 3: The spinning SC attitude estimation errors. 
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Spinner Rate Errors (deg/sec) 



Fig. 3: The spinning SC attitude estimation errors. 


V .2 The D-Filter Results 

The D-filter was run under the same conditions. The plots of the attitude and rate errors that are plotted on 
the same scale as the q-filter plots look similar to those shown in Figs. 3 and 4, respectively. The mean and 
standard deviation of the attitude and the angular rate errors are presented, respectively, in Table III and 
Table IV. 

Table III: Attitude error mean and standard deviation using the D-filter (deg). 



Component 

Mean error 

Standard Deviation 

ATTITUDE 

ERRORS 

X 

-0.094162 

0.036783 

y 

-0.169427 

0.018187 

z 

0.032968 

0.022750 


Table IV: Rate error mean and standard deviation using the D-filter (deg/sec). 



Component 

Mean error 

Standard Deviation 

RATE 

ERRORS 

X 

0.000006 

0.000004 

y 

0.000014 

0.000015 

Z 

- 0.000000 

0.000000 
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VI. CONCLUSIONS 


This paper treats the problem of recursive attitude and rate estimation of spinning spacecraft from vector 
measurement. The sensors considered are 3-axis magnetometers, sun sensors, star trackers, and horizon 
sensors. Error models suitable for Kalman filtering of these sensors were developed. Since usually spinning 
spacecraft carry no gyros, it was necessary to use Euler's equation of the angular dynamics of rigid bodies to 
estimate the angular rate. The use of the nonlinear Euler’s equation and the fact that quaternions are related 
to vector measurements in a nonlinear fashion turned the estimation problem into a nonlinear one. Following 
our good experience with the application of pseudo-linear Kalman filtering to similar problems, we adopted 
that filter in this work too. 

Two pseudo-linear Kalman filters were designed; namely, the q-filter and the D-filter. The former 
estimates attitude in terms of the attitude quaternion, and the latter does it in terms of the direction cosine 
matrix. Since usually vector measurements are not acquired simultaneously we used vector projection to 
project forward past measurements to coincide with present ones and thereby increase the robustness of the 
filters. 

A simpler Markov model is suggested as a replacement for Euler’s equation in the case where the vector 
measurements are obtained at high rates relative to the spacecraft angular rate. The performance of the two 
filters was examined using simulated data which confirmed the efficiency of the recursive algorithms in 
estimating attitude and rate of spinning spacecraft. 
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Appendix A: Sun Sensor Measurement Equation 2 

The geometry involved in the sun sensor measurement is presented in Fig. A. I 2 . 

Z b 



The known angles are the following: 

Measured : 

♦ 

Q 

9 

From Ephemeris : 

'F - arc length between the sun and earth. 

The computed unknown angles are: 

O - dihedral angle between the sun line and the center of earth. 

T| - Nadir Angle, the arc angle between the Z b and the direction to the center of earth. 

From Fig. (A. 1) it is seen that the computation of ® i s rather trivial because 


- dihedral angle between the sun line and earth horizontal crossing. 

- earth width (dihedral angle from earth-in to earth-out). 

- arc length between the sun line and Z b . 
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Q 

O = (J)+ — 

2 

The computation of T| is more involved. We start by using the following connection [2] 

cos W = cos <p • cos r\ + sin (p ■ sin r| • cos O 
In order to solve the last equation for r\ define 

a = cos(cp) 
b = sin(cp)cos(0) 
c = cosQ¥) 


then Eq. (A.2) becomes 


c = a cos(ti) + bsin(r|) 


Divide this equation by Va 2 +b 2 to obtain 


c a b 

r cos(r|) + y sin(r|) 


) / 

Va 2 +b 2 Va 2 +b 2 Va 2 +1 


Define an angle X such that 


and 


sin(X) = 


V a 2 + b 2 


cos(A.) = 

The defmition of this angle is explained in Fig. (A.2). 


Va 2 + b 2 



Fig. A.2: Definition of the angle 


l as 


Also let A = c/>/a 2 -fb 2 then Eq. (A.4) can be written ; 

A = sin(A.) cos(t]) + cos(X) sin(ri) 

Using the formula for the sine of the sum of two angles, the last equation can be written as 

A = sin(A. + rj) 
r| = arcsin(A) - X 


therefore 

which can be written as 
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(A.l) 

(A.2) 

(A. 3) 
(A.4) 

(A.5) 

(A. 6.1) 
(A. 6.2) 


(A.7) 

(A.8) 

(A.9) 



fa> 


r| = arcsm(A) - arctan — 

\bj 

= arcsin 


Substituting the values for a, b and c given in Eq. (A.3) yields 


■b 2 J 


- arctan I 


Va + 
the final solution for 


b 

n 


r| = arcsm 


cosY 


^ a/cos 2 <p + sin 2 (p- cos 2 <£> J 


-arctan - 


coscp 


V sin cp ■ cos O y 


and using Eq. (A. 1) the last equation becomes 
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r| = arcsm 


cos'? 


^/cos 2 (p + sin 2 (p-cos 2 ((() + Q/2) 


- arctan 


COS(p 


sin(p-cos(<|> + Q/2) y 


With T| on hand we can compute the unit nadir vector in boddy coordinates as follows 


1 


hb 


sin r| • cos(<j> + Q / 2 + 9) 
sin r| • sin(<J) + Q / 2 + 9) 
cosr| 


(A. 10) 


(A. 11) 


(A. 12) 


(A. 13) 
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